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Quantum Monte Carlo (QMC) techniques are used to calculate the one-body density matrix and excitation 
energies for the valence electrons of bulk silicon. The one-body density matrix and energies are obtained from a 
Slater- Jastrow wave function with a determinant of local density approximation (LDA) orbitals. The QMC density 
matrix evaluated in a basis of LDA orbitals is strongly diagonally dominant. The natural orbitals obtained by 
diagonalizing the QMC density matrix resemble the LDA orbitals very closely. Replacing the determinant of 
LDA orbitals in the wave function by a determinant of natural orbitals makes no significant difference to the 
quality of the wave function's nodal surface, leaving the diffusion Monte Carlo energy unchanged. The Extended 
Koopmans' Theorem for correlated wave functions is used to calculate excitation energies for silicon, which are in 
reasonable agreement with the available experimental data. A diagonal approximation to the theorem, evaluated 
in the basis of LDA orbitals, works quite well for both the quasihole and quasielectron states. We have found 
that this approximation has an advantageous scaling with system size, allowing more efficient studies of larger 
systems. 

PACS: 71.10.-w, 71.20.-b, 71.55.Cn 
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I. INTRODUCTION 

The two most common, practical quantum Monte 
Carlo (QMC) methods for realistic systems are the vari- 
ational quantum Monte Carlo r£VMC)t!rEl and diffusion 
quantum Monte Carlo (DMC)BB methods. In VMC, 
expectation values are computed with an approximate 
many-body trial wave function. In DMC, imaginary time 
evolution of the many-body Schrodingcr equation in prin- 
ciple gives exact results, although in practice one needs to 
make the "fixed-node approximation" to account for the 
antisymmetry of the many-electron wave function. In the 
fixed node approximation, the nodes of the propagated 
wave function are restricted to those of the trial wave 
function. The accuracy of this approximation is central 
to DMC simulations of many-electron systems. One of 
the aims of our work is to investigate the effectiveness of 
this approximation for extended systems, with the long 
term goal of obtaining better trial wave functions. 

In this paper we calculate the one-body density ma- 
trix for the valence electrons of silicon within the VMC 
framework, and obtain the natural orbitals which diago- 
nalize the density matrix. These calculations require the 
whole of the density matrix throughout all of the six- 
dimensional space r x r', not just at a few points in space 
as has been obtained before. To our knowledge this is the 
first time that the one-body density matrix and natural 
orbitals have been obtained for an extended, inhomo- 
geneous, interacting electron system. Recent evidence 
has suggested^ that a determinant of natural orbitals 
may give a better nodal surface than a determinant of 
Hartree-Fock (HF) orbitals. Our results show that a de- 
terminant of natural orbitals has a similar quality nodal 
surface to a determinant of LDA orbitals for bulk silicon. 
In a separate calculation we find that a determinant of 



LDA orbitals has a slightly better nodal surface than a 
determinant of HF orbitals. 

There is considerable interest in calculating excitation 
energies using QMC techniques. Excitation energies may 
be obtained by analyzing DMC decay curves,Ertl but this 
method has not proven very useful due to the large sta- 
tistical noise. Furthermore, as the quality of the ground 
state trial wave function improves, less information about 
excited states is obtained. A combination of ground and 
excited state wave functions must then be used to obtain 
upper bounds for the excitation energies. Direct methods 
for calculating excitation energies have met with more 
success. Mitas and Martin have calculated an excitation 
energy in a molecular nitrogen solid by performing DMC 
calculations for the ground and excited states.Q Mitas has 
also reported similar calculations for two excitation ener- 
gies in diamonds Recentlyu we used the same method to 
calculate 27 excitation energies in silicon, obtaining very 
good agreement with experiment for the low lying excita- 
tion energies, while the energies of the higher lying exci- 
tations were somewhat too large. In this paper we calcu- 
late excitation energies using a different approach. Jfece 
we use the "Extended Koopmans' Theorem" (EKT) Ji3lLil 
which derives from quantum chemistry, and involves the 
one-body density matrix. We have applied this theorem 
within VMC to calculate the excitation energies of silicon 
at four inequivalent k-points within the Brillouin zone. 
The energies are in good agreement with the available 
experimental data with a level of agreement similar to 
direct excitation calculations. We also test the diago- 
nal approximation to the EKT evaluated using the LDA 
orbitals, which was used previousi¥_.to estimate quasi- 
hole energies in siliconllj and NiO.Lj We find that the 
approximation performs well in silicon and that it has an 
advantageous scaling with system size. This allows more 
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efficient studies of excitations in large systems than are 
possible with existing direct techniques. 

The layout of this paper is as follows. In section n] we 
briefly describe the QMC techniques used in our calcula- 
tions, including the Hamiltonian, the trial wave function 
and the relevance of natural orbitals to QMC calcula- 
tions. In section III we present and discuss our results 



for the one-body de nsity matrix and the natural orbitals 
of silicon. In section IV we describe the Extended Koop- 
mans' Theorem and its application to the band structure 
of silicon. 



II. QMC SIMULATIONS OF SILICON 

In this section we briefly describe our QMC calcula- 
tions. For a more detailed discussiop-,oLthe methods we 
refer the reader to the literatureoF 11 ^ 



A. The Hamiltonian 

For this study we used an fee simulation cell, with pe- 
riodic boundary conditions, containing 54 Si 4+ ions and 
216 electrons. The Hamiltonian for our system, within 
the Born-Oppenheimer approximation, is 

i i a 



The positions of the N electrons in the supercell are 
denoted by and the ion locations are denoted by 
d Q . The electron-ion potential, v a , is modeled by a 
norm-conserving non-local pseudopotentialEEl obtained 
from atomic calculations performed within the local den- 
sity approximation (LDA) to density functional theory. 
The standard method for including the inter-particle 
Coulomb interactions in periodic systems is to use the 
Ewald interaction potential. We have found that this in- 
teraction gives rise to significant finite size errors, espe- 
cially for small simulation cells. Recently we introduced 
a new formulation of the electron-electron interaction for 
simulations using periodic boundary conditions which 
eliminates this problemO (hereafter referred to as the 
"cutoff interaction" ) . This interaction satisfies the con- 
ditions that (i) it gives the correct Hartree energy and (ii) 
it has the proper 1/r form for the interaction of an elec- 
tron with its exchange-correlation hole. (The Ewald in- 
teraction violates condition (ii).) Here we present results 
for excitation energies calculated with both the Ewald 
and cutoff interactions, using a wave function which was 
optimized using the cutoff interaction. For consistency 
one should use the same form of interaction between all 
the particles, but it turns out that if we apply our new 
interaction to a system of quantum mechanical electrons 



and classical ions then it reduces to using the Ewald in- 
teraction for the terms involving the ions while the cutoff 
interaction applies only to the electron-electron interac- 
tions. Note that the cutoff interaction is formulated in- 
dependently of QMC itself and jaiay be used with other 
techniques for periodic systems.ED 



B. The trial wave function 

The choice of trial wave function is of critical impor- 
tance for VMC and DMC calculations. We have used a 
standard Slater- Jastrow form: 

* T (ri, ...,rjv) = D T (ri, . . . , r« )D L (ri± +1 , ...,r N ) 

xexp(££iX(rO-E^«fo)) , (2) 

where the spin-up and spin-down Slater determinants, 
£> T and D x , are multiplied by a Jastrow factor which 
contains a one-body x function and two-body correla- 
tion factor, u. Our x-function has the full symmetry 
of the diamond structure and is expressed as a Fourier 
series containing 6 inequivalent, non-zero, parameters. 
We used spherically-, symmetric parallel and antiparal- 
lel spin u-functktnsjl3 which satisfy the electron-electron 
cusp conditionsEEl and contain a total of 16 parameters. 
The optimized parameter values w£re obtained by mini- 
mizing the variance of the energyJi3 

The spin-up and spin-down Slater determinants were 
formed from single-particle orbitals obtained from an 
LDA calculation employing the same pseudopotential as 
in the QMC calculations. The LDA orbitals were calcu- 
lated at the T-point of the simulation cell Brillouin zone 
using a plane wave basis set with an energy cutoff of 15 
Ry. Although the T-point scheme does not give opti- 
mal Brillouin zone samplingjH3 it does preserve the full 
symmetry of the system and allows comparison with a 
wider number of established results. The T-point of the 
simulation cell Brillouin zone unfolds to four inequiva- 
lent k-points in the primitive Brillouin zone. These are: 
(0,0,0) (the T-point), (0,0,|)^ (a point along the A axis, 
hereafter referred to as the A-point), (0,|,|) — (a point 
along the E axis, hereafter referred to as the S-point), 
and (!>!>|)^f (a point along the A axis, hereafter re- 
ferred to as the A-point). 

It is highly desirable to improve the quality of the trial 
wave functions used in QMC calculations. Improvements 
to trial wave functions can be classified into three types: 
(i) improvement of the Jastrow factor, (ii) using a lin- 
ear combination of determinants, and (iii) improvements 
in the orbitals forming the determinants. In this pa- 
per we will investigate a possible improvement of type 
(iii) , namely the use of natural orbitals. The question of 
which single-particle orbitals lead to the best approxima- 
tion to the exact many-body wave function is still open. 
Furthermore, this choice fixes the nodal surface of the 
trial wave function and therefore determines the accuracy 
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of the fixed-node approximation. LDA and HF orbitals 
have been-used successfully in a number of atomic JiSEJ 
moleculaiE3n3 and solido QMC calculations, but so far it 
has not proved possible to perform a direct optimization 
of the single-particle orbitals of an extended system. A 
study of first row atoms and moleculesc3c3 showed that 
lower energies can be obtained in both VMC and DMC 
using a trial wave function containing several determi- 
nants obtained from a multi-configuration self-consistent 
field (MCSCF) calculation. However, a similar study 
for small silicon clusters found that trial wave functions 
containing a single determinant of natural orbitals com- 
puted within an MCSCF scheme gave better DMC re- 
sults than some multi-determinant wave functions u This 
result strongly suggests that the natural orbitals result in 
improved nodal surfaces, and motivates our calculation 
of the natural orbitals for bulk silicon. 

An expansion of a wave function in Slater determi- 
nants of natural orbitals requires a smaller number of 
terms fa r . a g iven accuracy than expansions using other 
orbitals.E3E3 Calculation of the natural orbitals is, how- 
ever, costly, . and , less expensive schemes such as natural 
pair orbitalsc3i23 have been proposed to improve con- 
vergence in quantum chemical calculations. It is not 
clear that orbitals arising in schemes designed to accel- 
erate convergence of configuration interaction (CI) cal- 
culations should give smaller fixed-node errors in DMC 
calculations than LDA or HF orbitals. However, as men- 
tioned above, there is some evidence to suggest that nat- 
ural orbitals have this property. Natural orbitals have 
not frequently been computed within fermion QMC, al- 
though VMC and DMC calculations of natural orbitals 
have been repeated for the ground states of the Li, C, 
and Ne atoms.EII No calculations of natural orbitals for 
realistic extended fermion systems have appeared in the 
literature to date, although for homogeneous systems the 
translational symmetry requires the natural orbitals to 
be plane waves. 

Systematic studies of multi-determinant wave func- 
tions in QMC are lacking for solids. It seems reason- 
able to assume that multi-determinant wave functions 
will have improved nodes, and therefore give a better 
representation of the exact wave function, but there is lit- 
tle direct evidence to support this. Multi-configurational 
approaches include correlation effects, but do so rela- 
tively inefficiently - large numbers of terms (configura- 
tions) are usually required to obtain a significant propor- 
tion of the correlation energy. This form of wave function 
is unattractive for QMC as we require an accurate rep- 
resentation of the wave function which can be rapidly 
evaluated. Therefore, we obtain the one-body density 
matrix and hence the natural orbitals from a VMC cal- 
culation using a correlated trial wave function, bypassing 
the need to determine them using a multi-determinantal 
calculation. 



C. VMC and DMC calculations 

In VMC we compute the expectation value of the 
Hamiltonian, H, or other operator, with a trial wave 
function, ^>t- This method gives a rigorous upper bound 
to the exact ground state energy. The Metropolis al- 
gorithm is used to generate electron configurations, R, 
distributed according to | v I't(R-)| 2 , and the energy cal- 
culation is performed by averaging the local energy, 
over this distribution. 

In our DMC calculations we use the short-time approx- 
imation for the Green's function with a time step of 0.015 
a.u., which has_been shown to give a small time-step er- 
ror in silicon.E3 Importance sampling is introduced via 
the trial wave function, ^t- We make the fixed node ap- 
proximation, restricting the nodes of the DMC solution 
to be those of the trial wave function. Approximately 
15xl0 3 statistically independent electron configurations 
were used and the acceptance/rejection ratio was greater 
than 99.9%. The computational cost of this method 
scales with the third power of the system size. Exact 
fermion techniques, such as the release node QMC and 
CI methods, have computational requirements increasing 
exponentially with the system size and are impractical for 
the system sizes used here. 

III. CALCULATION OF THE DENSITY MATRIX 
AND NATURAL ORBITALS 

A. Density matrix 

The one-body density matrixEl for a normalized wave 
function, ip, is defined as 

p(r,r')=N J r(r,r 2 ,-..,r N ) 

Xi/j(r',r 2 ,...,r N )dr 2 ...dr N . (3) 

To facilitate calculation we expand the density matrix in 
a basis of orbitals, 4>i, leading to 

p(r,r>)=Y,PiM*)<t>* j ( r ')- ( 4 ) 

We refer to the diagonal elements, pu, as the orbital oc- 
cupation numbers. For wave functions consisting of a 
single determinant, such as HF or LDA wave functions, 
the density matrix is idempotent (p = p 2 ) and takes the 
form of a sum over the occupied orbitals, i.e., 

N/2 

p(r,r')=2^(rK( r '), (5) 

so that the occupation numbers are 2 (including spin de- 
generacy) for occupied orbitals and for unoccupied or- 
bitals. 
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We write the matrix elements of the interacting density 
matrix, pij , as expectation values over the distribution 

M 2 : 



x \ip(n , . . . , r N ) | 2 dr'dri . . . dr N . 



(6) 



The permutation symmetry allows us to rewrite this in 
a way which is efficient for Monte Carlo evaluation. De- 
noting the average over the distribution, \ip\ 2 , as (. . -^p, 
the Monte Carlo expectation value is written as 




^(r n )^-(r')^ ! 



dr' 



(7) 



so that N values are accumulated at each step along the 
VMC walk. The integral over dr' is performed by sum- 
ming over a grid of uniform spacing whose origin is chosen 
randomly for each electron configuration. The same grid 
in r' is used for each term in Eq. ffl, which further reduces 
the computational cost. We tested a series of grid sizes 
for the r' integral, using identical configurations for each 
grid size to obtain a correlated sampling estimate of the 
difference between the integrals. We found that a grid 
containing 125 points in the simulation cell sampled the 
integral with sufficient accuracy. 

Provided that the density of points in the r' integral 
is kept constant, the statistical error in the individual 
elements of the density matrix for a given number of sta- 
tistically independent configurations is approximately in- 
dependent of system size. 

We used a basis set consisting of the lowest energy LDA 
orbitals at the 27 k-points in the primitive Brillouin zone. 
We tested the effect of varying the number of orbitals in 
the basis. We found that approximately 40 orbitals per k- 
point were sufficient, although to retain the symmetry we 
included all members of a degeneracy, so that the actual 
number used was either 39 or 40 orbitals, depending on 
the k-point. The normalization used in Eq. (3|) requires 
that 



N = Trp , 



(8) 



which provides a practical test for the completeness of 
the basis set. The total occupation of the matrix (Trp) 
was 215.9(2), which is within the statistical error of the 
number of electrons in the system, indicating complete- 
ness of the basis at the level of the statistical accuracy 
obtained. 

For a given number of Monte Carlo moves, the best 
statistics are obtained by accumulating all non-zero ma- 
trix elements and applying the symmetry afterwards. 
However it is computationally very expensive to accu- 
mulate all of them, and we found that a more efficient 
procedure was to accumulate only the independent non- 
zero matrix elements. The basis set of LDA orbitals are 



basis functions of the unitary irreducible representations 
of the symmetry group Oj l . Using .the "orthogonality 
condition for matrix representations"EJ we inferred that 
elements involving products of orbitals from inequivalent 
k-points and of differing representations are zero. We 
ensured that every occurrence of a given representation 
was identical, so that products between functions belong- 
ing to different rows were orthogonal. This procedure 
reduced the total number of independent and non-zero 
matrix elements, pij , from 42094 to 582 elements. 

These matrix elements were sampled using approxi- 
mately 6.6xl0 5 statistically independent configurations. 
The correlation lengths along the VMC walks of both the 
local energy and density matrix were found to be essen- 
tially the same. 



B. Results for the density matrix 

We found the matrix p^ to be very nearly diagonal, 
with little coupling between LDA orbitals. Double occu- 
pancy (spin-up and down) of orbitals is denoted by the 
value 2.0. The maximum difference between the inter- 
acting occupation number, pa, and the LDA occupation 
number was 0.0625(5), which occurred at the T25' state 
at the top of the valence band. The magnitude of the 
largest off-diagonal matrix element was 0.014(1), which 
is of similar order to the occupation number of the lowest 
unoccupied orbitals. The fractional errors in occupation 
numbers for orbitals of low occupation were large in com- 
parison to those of high occupation. We found that 97.6% 
of the total occupation of the density matrix is contained 
within the four occupied LDA bands at each k-point, and 
99.0% is obtained within the first ten bands. In Fig. |l| 
we plot the occupation numbers against the LDA band 
energies. The occupation numbers decrease almost lin- 
early with increasing LDA energy for both the occupied 
and unoccupied bands. 

In Fig. H we show the density matrix, p(r,r'), in the 
(110) plane, and the differences between the VMC and 
LDA matrices. The coordinate r is fixed at the center of 
a covalent bond, and r' ranges over the (110) plane pass- 
ing through the atomic positions. The density matrix 
consists of an asymmetric central peak, reduced in width 
along the bonding direction. A longer ranged structure 
is present in areas of high valence charge density, smaller 
by approximately one order of magnitude than the peak. 

The VMC value for the peak in the density matrix 
on the bond center at r = r' is 1.7% smaller than the 
LDA value. The VMC density matrix has a larger mag- 
nitude around the neighboring silicon ions than the LDA 
density matrix, which consequently has a slightly smaller 
range. We also examined the density matrix in intersti- 
tial regions, where we found more structure to be present. 
Again, the LDA and VMC results were very similar, with 
small differences between the two cases arising principally 
from the differing charge densities. 
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To investigate the effect of using a finite size simula- 
tion cell we compared the LDA density matrix computed 
for 3x3x3 and 4x4x4 k-point meshes, corresponding to 
simulation cells containing 54 and 128 atoms respectively. 
We found that the central peak was largely unchanged 
and the longer ranged structure was in qualitative agree- 
ment. The central maximum in the density matrix in 
Fig. H is at the point r = r' and its magnitude is directly 
proportional to the valence charge density at that point, 
which differed by 4.9% between the two simulation cell 
sizes. We expect the finite size effects in the QMC cal- 
culations broadly to follow those in the LDA, as we have 
found for the total energiesO 

In exact Kohn-Sham density functional theory the to- 
tal energy can be written entirely in terms of the one- 
body density matrix of the Kohn-Sham orbitals, whereas 
in a fully interacting system both the one-body density 
matrix and the pair-correlation function are required. 
Results for the pair-correlation function from accurate 
correlated wave functions and LDA calculations are ex- 
tremely different Exact Kohn-Sham density functional 
theory reproduces the exact charge density and therefore 
exactly reproduces the diagonal r = r' part of the density 
matrix. The off-diagonal part of the exact Kohn-Sham 
and interacting density matrices are not required to be 
the same. In silicon we expect the LDA to give a good 
approximation to the exact Kohn-Sham density matrix. 
For this system our results show that the entire density 
matrices are very similar in VMC and LDA. 

C. Natural Orbitals 

The natural orbitals were obtained by diagonalizing 
the density matrix in the basis of the LDA orbitals. An 
assessment of the statistical errors in the eigenvalues and 
eigenvectors was made by subjecting the matrix to ran- 
dom perturbations of order the statistical error. The 
eigenvalues varied by up to ±0.0004 on application of 
the small perturbations. All the calculated eigenvalues, 
Xi, of the_density matrix lie in the range < Xi < 2, as is 
required.^ Identical results were obtained when elements 
within statistical error of zero were explicitly zeroed. The 
overlap of the space occupied by the LDA orbitals and 
the corresponding natural orbitals is measured by the ab- 
solute value of the determinant of the matrix of overlaps 
between these two sets of vectors. This gave a value of 
0.9948, indicating that the spaces spanned are almost the 
same. 

The eigenvalues of the density matrix for the "oc- 
cupied" natural orbitals were very slightly larger than 
the corresponding matrix elements pa (by about 0.001). 
Consequently, the eigenvalues of the "unoccupied" nat- 
ural orbitals were very slightly decreased, so that Trp is 
invariant. Therefore a plot of the eigenvalues of the den- 
sity matrix would be indistinguishable from Fig. [l], which 
shows the diagonal elements of the density matrix in the 



basis of LDA orbitals. 



D. DMC Calculations 

As well as the LDA and VMC calculations, we per- 
formed fixed node DMC calculations with trial wave func- 
tions of the form of Eq. using LDA and natural 
orbitals to form the determinants. Re-optimization of 
the Jastrow and x functions to improve sampling effi- 
ciency in the DMC calculation was found to be unnec- 
essary. The resulting energies were —107.59 eV (LDA), 
-107.69(1) eV (VMC with LDA orbitals), -107.71(1) eV 
(VMC with natural orbitals), -108.10(1) eV (DMC with 
LDA orbitals), and -108.09(1) eV (DMC with natural 
orbitals). The VMC wave function appears to show a 
very slight improvement with natural orbitals compared 
with LDA orbitals. However, to within statistical accu- 
racy, the DMC energies obtained with LDA and natural 
orbitals are the same. This indicates that the nodal sur- 
faces given by the LDA and natural orbitals are of the 
same quality. 

E. DMC Comparison of LDA and HF Orbitals 

In light of these results it is interesting to compare the 
quality of the nodal surfaces obtained with LDA and HF 
orbitals, which are both commonly used in the determi- 
nantal parts of trial wave functions for QMC calculations. 

We investigated this by performing DMC calculations 
in silicon with an fee simulation cell containing 16 atoms. 
The smaller simulation cell enabled a large number of in- 
dependent configurations to be obtained rapidly. Wave 
functions expanded in a basis of atom-centeijed Gaus- 
sians were obtained from the HF and DFT codeEj CRYS- 
TAL95. We took special care to ensure that the LDA 
and HF calculations were done in equivalent ways to try 
and eliminate any bias in the comparison. A basis set 
of four uncontracted sp functions and one d polarization 
function per pseudo-atom was optimized separately for 
each calculation. The quality of the basis set is high - to 
obtain the same energy within a plane wave calculation 
would require a basis set cutoff of 12.5 Ry. We used the 
same non-local LDA pseudopotential as in our other cal- 
culations. In both calculations we used the same u and x 
functions and performed DMC simulations with an aver- 
age population of 640 walkers, performing approximately 
6.7 x 10 5 walker moves. We obtained DMC total energies 
of -107.488(3) eV per atom and -107.464(3) eV per atom 
for the LDA and HF guiding wave functions respectively, 
using the Ewald interaction in the many-body Hamilto- 
nian. 

The walker energies were approximately normally dis- 
tributed. Using a conventional t-test, the 95% confi- 
dence interval on the difference in energies obtained was 
0.002 — 0.046 eV per atom, showing that for this system 
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it is very likely that the DMC energy from a determinant 
of LDA orbitals is lower than that from a determinant of 
HF orbitals. Therefore, for this system, a determinant 
of LDA orbitals has a marginally better nodal structure 
than a determinant of HF orbitals. 



IV. EXCITATION ENERGIES 

A. Excited State Calculations 

The calculation of excited state energies in solids using 
QMC methods is a fairly new area of research. Signifi- 
cant successes have been achieved using direct methods, 
in which separate QMC calculations are performed for 
the ground and excited states, and the excitation en- 
ergy is calculated as the energy difference ilu In these 
direct methods a QMC calculation must be performed 
for each excitation. In contrast, for the method described 
here a large number of excitation energies are obtained 
from a single QMC calculation involving averages over 
the ground state wave function. 



B. The Extended Koopmans' Theorem 

Our method for determining excitation energies cor- 
responds to a QMC formulation of the Extended Koop- 
mans' Theorem IEKT) derived independently by Smith, 
Day and GarrodEj and by Morrell, Parr and Levy.EiJ The. 
EKT is closely related to the earlier work of FeynmanEll 
on calculating excitation energies in the superfluid state 
of He 4 , although the quantum chemists appear to have 
developed the theory independently. The EKT has been 
shown to give verp-|good excitation energies for simple 
molecular systems J 3 .?! and has been applied to atomic and 
diatomic systems P 9 " 4 ^ It appears particularly well suited 
to QMC calculations in which explicitly correlated many- 
body wave functions are used. Here we review the deriva- 
tion following Rcf. nfl and present our QMC formulation. 



1. Valence Band Energies 

In this method the band energies are calculated as ion- 
ization energies. We start with an approximation to the 
normalized ground state wave function, ip . The wave 
function for the N-l electron system is approximated by 
the Ansatz of eliminating an orbital from ip N : 



^ iV - i (r 2 ,...,r Af ) = / u* v (r 1 )ip"(T 1 ,...,r N )dr 1 . (9) 

The valence orbital to be eliminated, u v , will be de- 
termined variationally. This Ansatz is reminiscent of 
a quasiparticle wave function for the excited state, al- 
though the formulation is for N-l particle eigenstates of 



the Hamiltonian. Expressing Eq. (^|) in second quantiza- 
tion yields 



[0 



N-l 



>= O v \ip N > 



(10) 



where O v is the destruction operator for the state u v . 

The ionization energy is given by the difference in the 
expectation values of the Hamiltonian pGalculated with 
the N and iV-1 electron wave functionsillilci 



N . {O v iJj n \H\O v iP n ) 



(O v i> N \O v i> N ) 



(11) 



If tp is an eigenfunction of H, Eq. (|Tl]) may be written 
as 



^ N \6l[H,d v ]\4, N ) 



(12) 



The denominator in Eq. ([12]) is the one-body density ma- 
trix. We now expand in a set of orbitals, {fa}, so that 
u v (r) — Ci V <pi(r), and O v — J2ci V a,i, where dj is the 
destruction operator for fa. The condition for a station- 
ary value of e v generates a secular equation 



(V - e v S v )c v = . 



(13) 



The matrix S v is the one-body density matrix, and the 
elements of V are = {ip N \a}-[H, di]\tp N ), where 



V%=Nl fa(r 1 )fa J (r , )r(ri,...,r N ) 

xHii/j(t' ',r 2 , . . . ,r N )dr'dri . . .dr N 



(14) 



Hi consists of the terms in the A^-electron Hamiltonian 
of Eq. (Q) involving coordinate ri , so that 



N 

Hi = hi +y^v(r 1 ,r j ) , 



(15) 



where h consists of the one-body kinetic energy opera- 
tor and ionic potential, including both local and non- 
local pseudopotential components, and v is the electron- 
electron interaction potential. 

If a HF wave function is used for ip N and the density 
matrix is expanded in a basis set of HF orbitals, then 
Vjj reduces to a matrix with the HF A^-particle eigen- 
values on the occupied part of the diagonal, and zeros 
everywhere else. The resulting excitation energies axe. 
those given by the well-known Koopmans' theorem.E£l 
The contents of the EKT method are now reasonably 
clear. The method consists of a quasiparticle-like Ansatz 
for the wave function of the N-l particle system, which 
is used to calculate the ionization energies of the system. 
Electron correlations are included, but no allowance is 
made for relaxation of the other orbitals in the presence 
of the excitation. Although this relaxation can be im- 
portant in small systems, it is expected to be much less 
important for excitations in extended systems such as the 
silicon crystal studied here. 
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2. Conduction Band Energies 

An analogous theory exists for the conduction band 
energies. The wave function for the N+l electron system 
is approximated by the Ansatz of adding an orbital to 



^ iV+i (r , ...,r N ) = Au c (roW" (r x , . . . ,r N ) , (16) 

where A is the antisymmetrizer and the orbital u c is to 
be determined variationally. In second quantization we 
have 

^n+i >= 6l\ip N > . (17) 
The excitation energies e c arc defined by 

{il> N \d c [H,6t}\i; N ) 



{>>p N \6M\i> N ) 



(18) 



Expanding in a set of orbitals gives u c (r) — X) Cjc^i ( r ) 
and 0\ = ^Ci C a\. The coefficients c; c are the solutions 
of the secular equation 



(V c - e c S c )c c = , 



(19) 



where the matrix elements are V£ = (ip N \&i[H, Oj]\ip N ) 

and Sf. = (^ja<a]|^> = % - Pij ■ 

It is instructive to introduce a new potential with ma- 
trix elements Vij = V% + Vfi: 

Vij = J (f)i(r )ho(t>*{ro)dr 
+ N J 0i(ro)^(ro)u(r o ,ri) 

x |<0(n, . . . , rAr)| 2 dr dri . . . dr N 
- A/0 i (r o )0*(ri)i;(ro,ri)V'(ro,r2, ...,r N ) 

Xip*(r 1 ,..,,r N )dr dr 1 ...dr N . (20) 

If we use a HF wave function and expand in a basis of HF 
orbitals, Vij reduces to a matrix with the HF AT-particle 
eigenvalues on the diagonal, and zeros everywhere else. 
In this case one can readily identify the second and third 
terms in Eq. (|2FJ) as, respectively, the Hartree and ex- 
change terms. (Similarly Vy reduces to the HF energy 
eigenvalues on the unoccupied part of the diagonal.) 



3. VMC formulation of the EKT 



We accumulate the matrix elements of and Vij, 
subsequently forming the matrix V$ — Vij — V^. The 
matrix elements, V^, are given by 



V2 = N MnWAr') 



, H^irx, ...,t n ) 



t(j(r 1 , ...,r N ) 



ip(r',r 2 , ...,r N ) 2 
x — h0(ri, . . . ,rjv) \ dr dn . ..dr N . (21) 



As before, we use the permutation symmetry to write 
this as 



V v - 



, Hnip^x, ...,r N ) 



i/)(ti, ...,r N ) 
tp(. . . ,r„, . . .) 



dr 



(22) 



so that N values arc accumulated at each step. The terms 
contributing to H n tp/tjj are already available in a VMC 
calculation, allowing V^ to be accumulated with virtually 
no additional cost beyond that required for the density 
matrix. An analogous VMC formulation for calculating 
Vij was used. The sing le p article terms, ho, appearing 
in the first term in Eq. (|2C|) are evaluated directly with- 
out using Monte Carlo integration. We found that using 
Monte Carlo integration for all the terms resulted in a 
small increase in the variance of the matrix elements. 
Therefore we prefer to calculate the hg terms directly. 
The matrix elements and Vij were accumulated at 
the same time as the elements pij . The full crystal sym- 
metry was again used, which reduces the number of ma- 
trix elements which must be accumulated and ensures 
the correct symmetry in the presence of statistical noise. 



C. Results for excitation energies 

1. Full EKT results 

The correlation lengths along the VMC walks for V^ 
and V,j were found to be similar to those of the local 
energy and density matrix, and the distribution of sta- 
tistical errors was similar to that for the density matrix. 
The elements of V^ and V£ with the largest statistical 
errors were the diagonal elements. The statistical error 
bars on these elements are estimated to be ±0.2 eV. 

The matrix Eqs. ^ and ^ were diagonalized using a 
generalized eigenvalue solver. The ratios of the statistical 
error bars to the mean values were significantly larger for 
the diagonal elements of V^ and V£ than for p^. Just 
as for the density matrix, small perturbations were ap- 
plied to estimate the statistical errors in the eigenvalues 
and eigenvectors. The numerical stability of the diag- 
onalization was improved by explicitly zeroing elements 
of VZ and V c 4 within statistical error of zero. The accu- 
racy of the eigenvalues was further verified by gradually 
increasing the number of bands in the diagonalization 
procedure. The valence and low lying conduction band 
energies were stable to within ±0.4 eV. 

The resulting band energies are given in Table | and 
m Fig. |, with the energy at the top of the valence band 
set to zero. Of the k-points computed here, the avail- 
able experimental data (Exp) is limited to the T-point. 
Because of Ais we also give empirical pseudopotential 
(Emp) dataCJ in Table | and in Fig. || which should pro- 
vide a good interpolation between this data. Results for 
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the Ewald and cutoff interactions are in good agreement, 
indicating that the finite size errors are not significant at 
the level of statistical accuracy obtained here. The ener- 
gies are in good qualitative agreement with the empirical 
data at all k-points except for the upper Ei valence band 
state, and the A2' conduction band state. The source of 
the error is principally the value of the diagonal matrix 
elements for these states, V~ and V~ (see next section). 

The EKT band energies at the T-point are in good 
agreement with the available experimental data, and are 
also in good agreement with the DMC data from direct 
calculations of the excitation energies The EKT valence 
band width of 12.9(6) eV is smaller than the value of 
13.6(3) eV obtained from the DMC calculations and is in 
good agreement with the experimental value of 12.5(6) 
eV. In comparison the HF data shows the well known 
overestimation of band gaps and band widths which is 
due to the neglect of correlation energy, the LD A gives ex- 
cellent valence band energies while the conduction bands 
are too low in energy by 0.7-1.0 eV, and the GW data is 
in very good agreement with experiment. 

The EKT is a formulation for the eigenstates of the 
jV-1 and N+l electron systems, while the direct method 
is aimed at calculating the eigenstates of the N electron 
system. The EKT and direct results should therefore dif- 
fer by the exciton binding energy, but this energy is small 
for silicon and cannot be resolved at the level of statisti- 
cal accuracy obtained. Depending on the application one 
would like to be able to choose whether to include exci- 
tonic effects, so that it is advantageous to have the dif- 
ferent QMC techniques available. Clearly further refine- 
ment of excited state QMC methods is required, but the 
results from the EKT and direct approaches are promis- 
ing for the study of more strongly correlated systems, for 
which the LDA and GW methods are less reliable. 



2. Diagonal approximations 

If we neglect the off-diagonal element of V£ , V£ , and 
Pij then the valence and conduction band energies can 
be approximated by 

DEKT _ v ii DEKT _ v ii / 9 o\ 
iv — : e ic — -I _ ' \ A0 ) 

Pit L pa 

where the superscript DEKT denotes the diagonal ap- 
proximation to the EKT. This approximation has been 
used within VMC by Fahy et alJiSpto calculate quasihole 
energies in silicon and by TanakalU to calculate quasi- 
hole energies in NiO. This approximation has the com- 
putational advantage that far fewer matrix elements are 
required, and also that the problems of statistical noise 
are reduced, because the values of only two matrix ele- 
ments enter the calculation of each band energy. How- 
ever the results are basis set dependent, and could differ 
significantly between, for example, a HF and LDA basis. 



If we use a HF wave function and expand in a basis set 
of HF orbitals then the matrices , , and p^ are all 
diagonal, and consequently the full EKT and the DEKT 
are equivalent. In general, for correlated wave functions, 
the DEKT gives neither an upper nor lower bound to the 
energy obtained from the full EKT. Comparison with the 
data from the DEKT is also made in Table |. The DEKT 
values are close to the full EKT values, including the two 
cases mentioned above where the agreement with experi- 
ment is poor. The DEKT works quite well for the valence 
bands and slightly less well for the conduction bands, be- 
cause the off-diagonal matrix elements of Vfi coupling the 
unoccupied states are more significant. In similar VMC 
calculations for silicon using the DEKT— and a simulation 
cell containing 64 electrons, Fahy et al.E3 calculated a va- 
lence band width of 14.5(4) eV, which is larger than both 
our value of 13.3(2) eV, and the experimental value of 
12.5(6) eV. Fahy et al.Ej also obtained occupation num- 
bers of 1.96(4) and 1.92(3) for the Ti and T25' valence 
band states, respectively, which are within error bars of 
our results of 1.9817(2) and 1.9375(5). 

The scaling with system size of the diagonal approxi- 
mation is very advantageous compared with direct evalu- 
ation of excitation energies. In direct methods, the frac- 
tional energy change due to the promotion of an electron 
is considered. This energy change is inversely propor- 
tional to the number of electrons in the system, i.e., a 
'■h' effect. The precision of the calculation must there- 
fore be sufficient to resolve the energy change amid the 
statistical noise. This requirement is very challenging 
for small band gap materials, as the system must also 
be sufficiently large to approximate the infinite solid. In 
the DEKT, the band energies are computed directly as 
averages over the square of the ground state wave func- 
tion, rather than as the difference between averages over 
the squares of the ground and excited state wave func- 
tions. This greatly improves the sampling statistics. We 
have found that the errors in V^, V£ and pa scale as the 
inverse of the square root of the number of independent 
electron configurations, independently of the system size. 
For sufficiently large systems, it therefore requires fewer 
statistically independent samples in the DEKT to obtain 
excitation energies to a given statistical accuracy than 
with direct methods. The scaling behavior of the full 
EKT is more difficult to analyze because we are required 
to diagonalize noisy matrix equations, where the number 
of off-diagonal matrix elements increases as the square of 
the system size, and where there are certainly statistical 
correlations between matrix elements. 



V. CONCLUSIONS 

We have calculated the one-body density matrix for the 
valence electrons of bulk silicon using QMC techniques. 

In real space, the VMC and LDA density matrices are 
very similar, the greatest differences being due to the 
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differing charge densities in each case. The natural or- 
bitals, obtained by diagonalizing the density matrix, very 
closely resemble the LDA orbitals. The occupation num- 
bers of the natural orbitals differ significantly from the 
non-interacting values, reducing linearly with increasing 
energy for the LDA occupied bands. The occupations are 
about 3% lower than the non-interacting value near the 
top of the valence band, and above the Fermi level, the 
occupation numbers fall slowly to zero. 

A DMC calculation for the ground state energy of sili- 
con using a trial wave function containing a determinant 
of natural orbitals gives an energy which is almost identi- 
cal to that obtained using a determinant of LDA orbitals. 
This shows that the quality of the nodal surfaces is al- 
most identical in each case. We used DMC calculations to 
compare the quality of the nodal surfaces obtained with 
LDA and HF orbitals, finding that the LDA orbitals gave 
a slightly lower DMC energy, indicating that for the sys- 
tem studied the nodal surface of a determinant of LDA 
orbitals is slightly better. We, note that a previous DMC 
study of small silicon clustersa found that natural orbitals 
obtained from a multi-configurational Hartree-Fock cal- 
culation gave a better nodal surface than a single deter- 
minant of HF orbitals. 

We have calculated excitation energies in silicon using 
an extension of Koopmans' theorem applicable to corre- 
lated wave functions. The Monte Carlo formulation is 
very similar to that required to obtain the density ma- 
trix. The resulting band energies are in good agreement 
with the available experimental data. 

The success of the VMC-EKT relies on a cancelation 
of errors between the ground and excited state energies. 
The wave functions for the excited states contain the vari- 
ational freedom of the orbitals u v and u c . This varia- 
tional freedom in the orbitals reduces the energy of the 
N-l and N+l electron states and improves the agree- 
ment with experiment. In the diagonal approximation to 
the EKT (DEKT), the u v and u c orbitals are fixed and 
there is no variational freedom in the excited state wave 
functions. We have found that the DEKT works quite 
well for silicon using LDA orbitals for u v and u c . The di- 
agonal approximation is exact within HF theory, so that 
we expect it to be a good approximation for weakly cor- 
related systems. 

Greater accuracy could be obtained with more accu- 
rate trial functions, or using DMC in the calculation of 
Vy, Vfj, and pij. In comparison with direct methods 
of calculating excitation energiesJ3i3 the EKT has the 
advantage that only a single calculation involving the 
ground state is required to obtain many excitation en- 
ergies. 

The EKT involves assumptions about the nature of 
the excited state wave functions, but nevertheless is a 
practical method for calculating excitation energies in- 
cluding correlation effects for relatively weakly correlated 
systems. The diagonal approximation to the EKT has 
a very advantageous scaling with system size compared 



with direct QMC calculations. This scaling thereby al- 
lows the study of excitation energies in larger systems 
with a greater efficiency than is possible with direct tech- 
niques. 

This work was supported by the Engineering and Phys- 
ical Sciences Research Council (UK). Our calculations are 
performed on the CRAY-T3D at the Edinburgh Parallel 
Computing Centre, and the Hitachi SR2201 located at 
the University of Cambridge High Performance Comput- 
ing Facility. 
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FIG. 1. Occupation numbers (diagonal elements of the den- 
sity matrix in the basis of LDA orbitals) plotted against the 
LDA energy for each k-point. These are the r-point(o), A 
(O), E (v)i an d A ( D )i points on the A, E, and A axes re- 
spectively. The statistical error bars are approximately equal 
to the sizes of the symbols for the conduction band states and 
are about 5 times smaller than the symbols for the valence 
band states. 
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FIG. 2. (a) The VMC one-body density matrix, 
pvMc(r, r') and (b) p V Mc(r,r') - Plda (r, r') , in the (110) 
plane passing through the atoms with r fixed at the bond 
center, p is normalized such that p(r, r) = n(r), the charge 
density at the bond center. The silicon atoms and bonds are 
shown schematically. 
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FIG. 3. The band structure of silicon, as computed within 

the VMC Extended Koopmans' Theorem using the cutoff in- 
teraction (□) and the Ewald interaction (o), and within the 
LDA (x). The results for the cutoff interaction are shown 
with statistical error bars. The statistical error bars for the 
Ewald results have been omitted for clarity, but they are the 
same size as for the cutoff interaction. As a guide to the eye, 
the empirical pseudopotential data of [ ^] is shown (solid 
lines), which is in good agreement with the available experi- 
mental data. 
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TABLE I. Band energies of silicon in eV. a,b- VMC-EKT energies using the cutoff and Ewald 
interactions, respectively. The statistical error bars are ±0.4 eV. c,d- Diagonal approximation to 
the VMC-EKT energies using the cutoff and Ewald interactions, respectively. The statistical error 
bars are ±0.2 eV. e- Direct DMC calculations, with statistical error bars of ±0.2 eV, from Ref. [j. 
f- Ref. g- This work, h- Ref. i- Ref. From the compilation given in Ref. H. 
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